function [err, data_set, base_targets, model, g_ind] = obj_func(x)
  
global bet0 bet1 tau a b w e  N_a N_w N_e w_max w_min alpha xi1 lb ub pi_e...
        sigmae elasts log_e_mu log_e_sigma skill_shock...
        base_rs shares wspread pe thetas pi_theta rho sig...
        unemp zeta  myslope nu xi eta  mu_a gamma b_init bet0_init bet1_init base_rs_init Pr share curve scale w_gamma

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Data Moments
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    % Housing moments: Share A 1980-2010,Share B 1980-2010, Rent Ratio A-B 1980, Rent Ratio B-C 1980, 
    Hs_moments_base = [0.1934594	0.21735	0.2278978	0.2508063	0.3006039	0.2494808	0.2294	0.2151092	1.2530   1.2769		];    
    Hs_moments = Hs_moments_base;
    
    % Gini 1980-2010
    gini_data = [0.3763; 0.4112; 0.4411; 0.4661];

    % Dissimilarity 1980-201
    dissim_data = [0.3341891; 0.3926; 0.4148; 0.4434701];

    % Education returns 1980-1990
    educ_return_data = [0.391, 0.549];

    % Education moments:
    % College share A,B,C and dissimilarity by education
    edu_data = [0.3348; 0.1791; 0.1131; 0.2495];
    
    % Rich share A, B, C
    rich_share_data = [0.4372;	0.2252;	0.0943];

    % Rank-rank correlation
    ch_rank_rank = 0.341;

    % Return to spillover 25th p
    ch_25 = 0.062;

    % Return to spillover 75th p
    ch_75 = 0.046;

    % Income 25th/75th p
    w_25_75 = 0.667;

    % Collect targets
    base_targets = [gini_data(1), dissim_data(1), ch_rank_rank, ch_25, ch_75, w_25_75];

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Model Parameters
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    b = x(1) ;
    
    xi =  abs(x(2)); 
    alpha = x(3);
    sig = abs(x(4));
    rho = abs(x(5));
    log_e_sigma = abs(x(6));
    gamma = x(14);
    tau = x(7);
    bet1 = (x(8));
    skill_shock = abs(x(9));
   
    % Preference Shocks
    theta_h = max(x(10), 1);
    theta_l = max(min(x(11), 1), 0); 

    % Preference Shocks Probility
    pi_theta = max(min(abs(x(12)), 1), 0);
    thetas = [theta_h, 1, theta_l];
    pi_theta = [1 - pi_theta, pi_theta];
     
    % Idiosyncratic shock variance
    sigmae = 1 /  abs(x(13))  ;  
    
    
    % These parameters can have arbitrary normalizations but aide
    % calibration
    base_rs = x(15); % Normalize rental rate in C so that other prices are relative to C.
    bet0 = (x(16));  %  Normalize bet0 so that maximum education is 1 in SS;
    scale = 1;  % Use this to scale the distribution of wages to 1980 median family income
    b = b*scale^(1+alpha);
    bet0 = bet0*scale^((1+alpha)/xi);
    bet1 = bet1*scale^((1+alpha)/xi);
    
    % Housing Parameters - Elasticities
    phi2 = x(end-1); 
    phi3 = x(end);
    phi1 = 1.75*3-phi2-phi3;

    
    % Elasticity
    elasts = [phi1, phi2, phi3 ]';
    Hs = ones(size(elasts));

    % Optional parameters
    pe = 0.0;    % Public education
    unemp = 0.0; % Unemployment benefits
    xi1 = 1;     % Curvature on spillover
    w_gamma = 1; % Curvature on warm-glow preferences
    
    % Housing Shares
    h_a = Hs_moments_base(1); 
    h_b = Hs_moments_base(5); 
    h_c = 1 - h_a - h_b;
    shares = [h_a; h_b; h_c];
       
    % Skill premium shock in SS
    eta = 1;

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Grids
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    N_at = 1;   % Grid size for temporary component of ability
    N_ar = 40;  % Grid size for AR(1) ability
    N_a = N_ar * N_at;
    log_ar_sigma = sqrt(sig^2 / (1 - rho^2)); % Implied std. deviation
    log_ar_mu = 0.0;
    log_at_mu = 0.0;
    log_at_sigma = 0.0;
    
    % Parameters governing random component of wages
    N_w = 125; % Grid size for Income
    log_e_mu = 0.0;
    log_w_mu = log(base_rs);
    log_w_sigma = 0.9;
    corr_wa = 0.0;
    N_e = 120; % Grid over mean zero wage innovations
    
    lb = 0.0; % Lower bound of education
    ub = Inf; % Upper bound of education
   
    % Define convergence
    tol = 1e-4;
    max_It = 150;

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Ability Process
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Grid and transition matrix for ability
    m_a = 4; % Std Deviations Away from SS
    [log_ar, Pr_ar] = tauchen(N_ar, log_ar_mu, rho, sig, m_a);
    ar = exp(log_ar);
    % Transitory component of ability
    at = 1.0; 
    pi_at = 1.0;

    % Combine into ability matrix
    a = kron(at, ar);
    a = a(:);
    log_a = log(a);
    log_a_mu = log_ar_mu + log_at_mu;
    % Std Deviation of log ability process
    log_a_sigma = sqrt(log_ar_sigma^2 + log_at_sigma^2);

    [~,a_ind] = sort(a);
    a = a(a_ind);

    Pr = repmat(kron(pi_at, Pr_ar), N_at, 1);
    Pr = Pr(a_ind, a_ind);
    norm_mean = dot((Pr'^50)*ones(N_a,1)/N_a, a);
      
    a = a / norm_mean ;
    base_rs = base_rs*scale;
    tau = tau*scale;

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Wage Grid
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    high_low = 4.5;
    w = zeros(N_w,1);
    w(end) = high_low * log_w_sigma;
    w(1) = log_w_sigma*-high_low;
    zstep = (w(end) - w(1)) / (N_w-1);
    w(2:end-1) = w(1) + zstep * (1:N_w-2);
    wspread = exp(w);
    w = w + log_w_mu;
    w = exp(w);
    w_min = min(w);
    w_max = max(w);
    log_w = log(w);
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Idiosyncratic shock of wages
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    e_min = logninv(0.00005, log_e_mu, log_e_sigma);
    e_max = logninv(0.99995, log_e_mu, log_e_sigma); 
    e = exp(log(e_min) + (log(e_max) - log(e_min)) * linspace(0, 1, N_e+1));
    e = e';
    pi_e = logncdf(e, log_e_mu, log_e_sigma);
    pi_mass = diff(pi_e); 
    pi_e = pi_mass / sum(pi_mass);
    e =  (e(1:end-1) + e(2:end))/2;
    e = e / dot(e, pi_e);
   
    
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Initial Disribution of wealth and ability
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    mu = [log_a_mu, log_w_mu];
    % Covariance Matrix
    Sigma = [log_a_sigma^2 corr_wa*log_w_sigma*log_a_sigma;
             corr_wa*log_w_sigma*log_a_sigma log_w_sigma^2];

    [m1,m2] = meshgrid(log_a, log_w);
    Prob = zeros(size(m1));

    for i = 1:size(m1, 1)
        for j = 1:size(m1, 2)
            Prob(i,j) = mvnpdf([m1(i,j), m2(i,j)],mu,Sigma);
        end
    end

    Prob = Prob / sum(Prob(:));   
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Initial guess for rental rate and spillover
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    Ss = [2.49569237238145;2.19773991908705;1.89032785122845]';
    Ss = Ss ./ dot(Ss, shares);
    rs = [base_rs*2, base_rs*1.5, base_rs]';
   
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Steady State
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    tic()
    [Probss, Sss, Rss] = steady_state(Prob, Ss.^xi1, Hs, rs, tol, max_It);
    toc()

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Summary Statistics
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Generate model moments
    data_set = generate_summary_statistics(Probss, Sss, Rss, Hs, tol, max_It);
    
    % Model Moments
    dissim = cell2mat(data_set(:,1));   % Dissimilarity     
    gini = cell2mat(data_set(:,2));     % Gini
    
    avg_inc = cell2mat(data_set(1,3));  % Average Income
    total_educ_return = cell2mat(data_set(:,8)); % Education Returns
    
    spillovers = cell2mat(data_set(:,9)); % Spillovers
    
    myrs = cell2mat(data_set(1,10)); % Rental rates
    
    col_share_model = cell2mat(data_set(1,16));     % College Shares
    dissim_ed_model = cell2mat(data_set(1,17));     % Dissimilarity education
    edu_model = [col_share_model dissim_ed_model];  % Education moments
    rich_share_model = cell2mat(data_set(1,14))';   % Rich Share
    
    % Compute some moments averaged between 1980-2000
    avg_T = 3; 
    population = 1.0893.^[0;1;2];   
    rankrank_avg = sum(population .* cell2mat(data_set(1:avg_T ,6)))/sum(population);   % Rank-rank
    spillover_25_avg =  sum(population .* spillovers(1:avg_T ,1))/sum(population);      % Return to spillover 25 p
    spillover_75_avg = sum(population .* spillovers(1:avg_T ,2))/sum(population);       % Return to spillover 75 p
    inc_ratio_avg = sum(population .* cell2mat(data_set(1:avg_T ,11)))/sum(population); % Income 25th/75th p
   
    
    % Identify the neighborhoods
    [~,idx] = sort(avg_inc);
    A_idx = idx(end);
    B_idx = idx(end-1);
    C_idx = idx(1);
    
    % Housing Shares
    Hs_all = cell2mat(data_set(:,7));
    A_series = Hs_all(1:4,idx(end))';
    B_series = Hs_all(1:4,idx(2))';
    C_series = 1 - A_series - B_series;
    
    % Housing moments
    hs_mom_model = [A_series, C_series, myrs(A_idx)/myrs(B_idx), myrs(B_idx)/myrs(C_idx)];
    Hs_moments(5:8) = 1 - Hs_moments(5:8) - Hs_moments(1:4);
    
    % Baseline model targets
    model = [gini(1), dissim(1), rankrank_avg, spillover_25_avg, spillover_75_avg, inc_ratio_avg];
  
    % Define weights
    % More important moments are upweighted by a factor 5
    weights = ones(size(base_targets));
    weights(1:2) = 5; 
    weights(3:5) = 5; 
    base_mom = abs(base_targets - model)./ base_targets  .* weights;
    
    % Educational returns moment;
    return_mom_weight = ones(1,2)*5;
    return_mom = abs(educ_return_data  - total_educ_return(1:2)') ./ educ_return_data  .* return_mom_weight;
    
    % Education moments;
    col_weights = ones(4,1);
    col_mom = abs(edu_model' - edu_data )./edu_data .* col_weights; 
    col_mom = col_mom(1:end-1); 
    
    % Rich share moments;    
    rich_share_weights = ones(3,1);
    rich_share_mom = abs(rich_share_model(:) - rich_share_data(:) )./rich_share_data(:) .* rich_share_weights; %
    
    % Housing moments
    HS_mom_weight = ones(length(hs_mom_model),1)*5; %5
    % (Last periods are not targeted, set weight to 0)
    HS_mom_weight(3:4)= 0;
    HS_mom_weight(7:8)= 0;
    neigh_mom = (abs(Hs_moments(:) - hs_mom_model(:)) )./ abs(Hs_moments(:)) .* HS_mom_weight; 
   
    % Collect all moments
    g_ind = [base_mom(:); neigh_mom; return_mom(:); col_mom; rich_share_mom(:)  ]; 
    
    % Objective to minimize
    err = dot(g_ind(:),g_ind(:)); 
    
    % Collect various model moments:
    model = [model, hs_mom_model, total_educ_return(1:2)', col_share_model, gini(end-1) , dissim(end-1)];
        
    disp(err)
    
end